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ABSTRACT 

We present the results from a series of collisionless iV-body simulations of major 
mergers of galaxy dark matter haloes with density profiles having either inner cusps 
or cores. Our simulations range from 2 x 10 5 to 10 7 particles, allowing us to probe the 
phase-space distribution of dark matter particles in the innermost regions (less than 
0.005 virial radii) of cold dark matter haloes, a subject of much recent debate. We 
find that a major merger of two cored haloes yields a cored halo and does not result 
in a cuspy profile seen in many cosmological simulations. This result is unchanged 
if we consider mergers with parent mass ratios of 3:1 instead of 1:1. Mergers of a 
cuspy halo with either a cored halo or a second cuspy halo of equal mass, on the 
other hand, produce cuspy haloes with a slightly reduced inner logarithmic slope. 
Cuspy haloes, once formed, therefore appear resilient to major mergers. We find the 
velocity structure of the remnants to be mildly anisotropic, with a Maxwellian velocity 
distribution near the centre but not in the outer portions of the final haloes. Violent 
relaxation is effective only during the early phase of mergers, with phase mixing likely 
to be the dominant relaxation process at late times. 
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1 INTRODUCTION 

In the currently favored cosmological constant plus cold dark 
matter (ACDM) model of cosmology, structure forms hier- 
archically: small dark matter haloes collapse at high red- 
shifts, while larger and more massive objects form through 
a series of minor mergers that accrete smaller mass haloes 
and major mergers with comparable mass haloes. Since dark 
matter haloes provide the g ravitational potentia l wells nec- 
essary for galaxy formation dWhite fc Reesll97Sft , a detailed 
understanding of the effects of mergers is necessary in order 
to test models of galaxy formation and evolution within the 
ACDM paradigm. 

Recent numerical studies indicate that the hierar- 
chical build-up process results in haloes with significant 
substructure populations: hundreds to thousands of sub- 
haloes to the smallest resolvable mass scales appear to re- 
side within galaxy sized haloes, comprising a pproximately 
10% of a halo's mass jChigna et all Il998t iKlypin et alJ 
119991: iMoore et al l Il999bl) . A number of recent papers 
have investigated the evolution and detailed properties 
of the substructure using both semi-analytic models and 
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numerical simulations. Depending on factors such as the 
relative mass, pericentric distance, and halo concentra- 
tion, the orbits and masses of the subhaloes can be 
strongly affect ed by dynamical fr i ction and tidal strip- 
ping ("e.g. [Velazquez fc White! Il999l: iTavlor fc Babull l200.lt 
lHavashi et al .120051: iTaffoni et al.ll2003l : lBenson et al]l2004Tl . 
Dense subhaloes that survive the tidal process can sink to 
the centre of the parent halo, impacting its inner stru cture 
(e.g. iDekel et alJl2003t iDekel. Devor. fc Hetzronll2003D . 

In contrast to minor mergers, where large haloes swal- 
low up smaller ones in a relatively gentle fashion, major 
mergers are more violent events involving systems of approx- 
imately equal mass that can be studied in detail with nu- 
merical simula tions only. The interest in major mergers has 
a long history. iToomre fc Toomrd dl972T) noted that many 
observed features of galaxy pairs can be explained by tidal 
interactions that could promote merging. They further spec- 
ulated that if such interactions were common enough, ellipti- 
cal galaxies might naturally be explained as the merger rem- 
nants of spiral galaxies. This suggestion has led to numerous 
studies of galaxy merging and has been quite successful in 
helping explain t he origin of elliptica ls. In some of the ear- 
liest simulations, White ( 1978], 119791) studied several types 
of mergers using 250 simulation particles per halo and con- 
cluded that products of mergers of cored galaxies were signif- 
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icantly more centrally concentrated than their progenitors, 
with more extended haloes having an approximate power- 
law density profiles p oc r~ 3 that resembled the de Vau- 
couleu rs' empirical profile fo r the surface-brightness of ellip- 
ticals lld e_Vauci3uleirrs]|l94d) ■ The higher resol ution simula- 
tions bv lPearce. Thomas, fc Couchmar] dl993l) with ~ 8000 
particles per halo also exhibited the r~ 3 behavior with evi- 
dence for homology. Simulations with multiple components 
representing dark matter haloes, stellar disks, and stellar 
bulges have been a source of much theoretical progress (e.g. 
iBarnes fc Hernauislll992l and references therein) and are es- 
sential for elucidating the dynamical coupling between the 
dark matter and stellar components in galaxies and for com- 
paring simulation results with observations. Pure dark halo 
simulations, however, are required for a controlled study of 
the non-baryonic dark matter that dominates the gravita- 
tional potential of seemingly all galactic systems. 

In this paper we perform high-resolution iV-body sim- 
ulations of mergers of two equal-mass dark matter haloes 
and study how the phase-space distribution of dark mat- 
ter particles is affected by such violent events and how the 
merger product relaxes into equilibrium. We pay particu- 
lar attention to the central regions of the haloes and ex- 
amine whether the initial inner structure of a dark mat- 
ter halo, either with a flat core or a sharp cusp, is dis- 
rupted or preserved by major collisions. Our study com- 
plements the numerous recent cosmological iV-body stud- 
ies of the properties and evolution of dark matter haloes, 
most of which have reached the somewhat surprising con- 
clusion that over a wide range of mass scales, these haloes 
are universally cuspy and have radial density profiles with 
p(r) oc r~ 7 (with 1.0 < 7 < 1.5) in the central region and 

oc 
1997 



_ the virial radius (Navarro. Frenk. fc White! 
IFukushiee fc Makinol Il997t iMoore et ail 



1999a 
20041 



Ghiena et alJ l2000l : IFukushiee. Kawai. fc Makinol 
The seemingly generic predict ion of an inner cusp agrees 
with earlier simulations (e.g. iDubinski fc Carlberd [1991) 
but is contrary to the inferred dark matter distribution 
in some dwarf and low surface brightness galaxies (e.g. 
iMcGaueh fc de BlokNl998t ISimon et al1l200A tho ugh also 
see|Swater^^tHnTl2OT3r as well as in spiral galaxies ISaluccil 
l200lT) and has lead to speculation that the ACDM model is 
in crisis. 

The simulations of major mergers performed in this pa- 
per provide a set of controlled experiments for us to address 
questions such as: How does the phase space density of dark 
matter haloes evolve and re-equilibrate during major merg- 
ers? Is the cuspy profile an attractor region of the phase 
space so that merger remnants of two non-cuspy (i.e. cored) 
haloes will end up cuspy? Do cuspy haloes, once formed (by 
any process), retain their inner structure after inevitable 
major mergers with other haloes? Unlike cosmological sim- 
ulations that are designed to incorporate many complicated 
processes at once, the simulations performed here allow us 
to separate out the effects of minor mergers and to focus on 
how major collisions change the spatial and velocity struc- 
tures of any initial haloes of our choice. Additionally, cos- 
mological parameters are unimportant for the simulations 
themselves since major mergers are local processes of two 
fully collapsed objects. In this case, all we need is the initial 
phase space distribution of each galaxy halo and the initial 
relative positions and velocities of the haloes themselves. By 



simulating collisions of two galaxy haloes with high spatial 
resolution and up to 5 million particles (for each halo), we 
are able to follow their phase-space evolution and the den- 
sity profiles of the remnants down to ~ 0.004 V2oo in our 
highest resolution simulation. 

Section 2 of this paper describes the initial setup and 
the parameters used in the numerical simulations. We also 
discuss in some detail the tests performed to quantify the 
numerical effects of two-body relaxation on the central re- 
gions of the haloes and the technical issue of how to deter- 
mine the centre of a merged halo. Section 3 contains our 
results from the ten production runs for the energy, density 
profiles, velocity distributions, and phase space structure of 
merger remnants. Section 4 discusses the energy distribution 
and the relaxation of a merged halo to an equilibrium state, 
with an emphasis on violent relaxation and phase mixing. 
Section 5 contains a summary of our results. 



2 PROCEDURE 

2.1 Initial Conditions 

We study two general types of inner density profiles for 
galaxy haloes: cored and cuspy. For the cored profile, we use 
an isot hermal sphere wi th a core radius and an exponential 
cut-off faernauisdli993h : 



p(r) 



aMexp(-r 2 /r c \ t ) 



(1) 



27r 3 / 2 r cut r 2 + r 2 OIC 
where M is the halo's mass, a is a normalization constant 
given in terms of x = r cu t/r COIC as a" 1 = 1 — y^rxe x erfc(x), 



- 1 = 1 

r cu t is the outer exponential fall-off radius and r, 



core is the 

core radius. For the cuspy haloes, we use the Navarro, Frenk , 
and White (NFW) profile jNavarro. Frenk. fc Whitelll997fl : 



p(r) = 



p c 8 



(r/r 3 )(l + r/r s 



(2) 



with scale radius r„ and characteristic density 8 given by 



5 = 



T200 _ 1 / 3M 2 Q0 X 
C ~~ C ^80071-pc 

200c 3 

3[ln(l + c) - c/(l + c) 



(3) 
(4) 



where c is the halo concentration parameter and r^oo is the 
radius at which the average interior density is 200 times the 
critical density of the universe p c = 3H 2 /8ttG. (We also 
adopt r2oo as a halo's virial radius and therefore M200 = 
8007rp c r-2oo/3 as its virial mass.) 

We gene rate our initial c onditions using the methods 
described bv iHernauisd l|l993h . The initial positions of the 
particles in a halo are drawn randomly from its density pro- 
file. All of our initial haloes are chosen to be 1O 12 M0. In the 
case of the NFW profile the total mass is formally infinite, 
so it is necessary to introduce a truncation radius for the 
halo. We choose to truncate our haloes at r = T2oo- For our 
choice of c = 10, the virial and truncation radius is 162.6 
h^ 1 kpc and the scale radius is r s = r2oo/c = 16.26 /i -1 kpc. 
For the cored profile, we choose r cor c = 8 h~ 1 kpc, r cut = 81 
h^ 1 kpc, and truncate the halo at 2r C ut. This choice of pa- 
rameters allows us to start with two types of galaxies, equal 
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in mass and radius, but with very different inner density 
profiles (core vs. cusp). 

For the initial velocities of particles in the frame of a 
given halo, we first use the Jeans equation (the first moment 
of the collisionless Boltzmann equation) to relate the density 
profile p(r) and potential <£>(r) to the velocity dispersions 
ilBinnev fe Tremaindll987l) : 



— h 2-pa r 

dr r 



where the velocity anisotropy is parameterized by 

„2 | 2 



(5) 



(6) 



and a\ is the dispersion of velocity component i. If we fur- 
ther assume initial velocity isotropy, j3(r) = 0, we can inte- 
grate Eqn. @ to yield the dispersion of the radial velocity 
as a function of radius, 



p(r) 



dRp(R)^. 



(T) 



From there, we randomly draw individual velocity com- 
ponents f rom a Gau s sian d istribution with variance <r 2 (r). 
Following iHernauistl il993fl . we impose the additional con- 
straint v < 0.95u csc , which insures that all particles are ini- 
tially bound to the halo. 

We set up the two colliding haloes to be just touching at 
the truncation radii initially. The initial relative speed of the 
two haloes is taken to be the speed one of the haloes would 
have if the two were in circ ular orbi t about their common 
centre of mass, i.e., v TC \ = 2yJ Gp>/r Te \, where p. is the reduced 
mass and r rc i is the distance between the halo centres. If the 
two haloes were point particles, this configuration would give 
^rci = v csc /\/2\ all our simulations are therefore relatively 
high-speed encounters for bound orbits. In this set-up, and 
with our choice of halo parameters, the two haloes first cross 
centres of mass at approximately 1.2 Gyr, then again at 2.2 
Gyr. This leaves ample time after the first encounter for the 
haloes to interact and evolve within the 5.0 Gyr run. 

Most of our simulations are head-on collisions with an 
impact parameter 6 = 0. We have also performed two sim- 
ulations in which the two halo centres are initially offset by 
b = 0.5r2oo in order to test the robustness of our conclusions. 
In terms of the dimensionless parameter 



A 



jm 1 ' 2 



(8) 



GM 5 /2 

for the total mass M, energy E, and orbital angular mo- 
mentum J, this corresponds to A = 0.042. In terms of the 
parameter e, defined as the ratio of the angular momentum 
to that needed for circular orbit at the same velocity and 
separation, this set-up yields e — 0.24. As a final test of the 
generality of our findings, we performed two simulations of 
collisions involving unequal mass haloes with cores. 

2.2 Simulation Parameters and Numerical Effects 

All of our simulations are done using GAD GET, a publicly 
availa ble N-body tree code iSpringel. Yoshida. fc White! 
l200lT) . GADGET is available in serial and massively parallel 
forms; we have tested both versions on a local linux clus- 
ter and the parallel version on the IBM SP2 computers at 



Table 1. Properties of ten simulations of halo mergers. All sim- 
ulations involve two 10 12 Mq haloes (except where noted), so 
M 20 o = 2 x 10 12 (7V 2O o /N p )M e 



Name 


b 1 

(rwo) 


N p 2 


N200 3 




r(t r ) 4 
(h- 1 kpc) 


Corel 







2 x 10 5 


1.58 x 


10° 


4.035 


Core2 





2 x 10 6 


1.58 x 


10° 


1.850 


Core3 


0.5 


2 x 10 6 


1.52 x 


10 6 


1.886 


Core4 5 





1.33 x 10 6 


1.16 x 


10° 


1.920 


Core5 5 





1.33 x 10 6 


1.09 x 


10° 


2.025 


Cuspl 





2 x 10 5 


1.34 x 


10 s 


3.595 


Cusp2 





2 x 10 6 


1.35 x 


10 6 


1.436 


Cusp3 





1 x 10 7 


6.72 x 


10° 


0.765 


Cusp4 


0.5 


2 x 10 6 


1.31 x 


10 6 


1.549 


Mixed 





2 x 10 6 


1.48 x 


10 s 


1.488 



1 impact parameter, in units of 7*200 f° r an initial halo 
2 total number of initial particles 

3 total number of particles within virial radius of remnant 
4 radius above which the two-body relaxation time exceeds 5 Gyr 
5 unequal mass (Mi = 3M 2 ) collision; Mi = 10 12 Af<r 



the National Energy Research Scientific Computing Cen- 
ter (NERSC). In these test runs, we use either 10 5 or 10 6 
particles and run for 5.0 Gyr. An important objective of 
the test runs is to check the level of stability against two- 
body relaxation for an isolated halo initially in equilibrium. 
Two-body relaxation is a numerical artifact of running sim- 
ulations where each simulation particle is on the order of 
10 60 - 10 75 times more massive than a typical elementary 
particle dark matter candidate. This leads to significant par- 
ticle encounters that one does not expect to find in a real- 
istic dark matter halo. As a result, the two-body relaxation 
timescale effectively puts a limit on the region of the halo 
that faithfully represents the phase space distribution in the 
collisionless limit. 

Though there has been some disagreement about 
what regions can be considered reliab le beyond the 
level of two-body relaxation (e.g. [M oore et al| |l99g|: 
[ Klvpin et alJl200lt iFukushige fc Makindl200ll : iPower et all 
120031: iDiemand et al]l2003l) . a conservative criterion is that 
the local two-body relaxation timescale t r , which we de- 
fine in terms of the circular velocity V c ii C = Vd rc (r), period 
T(r) = 2-Kr/Vch-c, and number of particles N (or equivalently 
mass M) interior to a radius r as 



t r (r) 



T(r 200) 



7T N 



N 



Kir. 



8 In N r 2 00 Kirc(»"200) 



4 In N V GM(r) ' 



(9) 



be at least as large a s the age of th e universe 
l|Fukushige fc Makind [20011: IPower et al.ll2003h . Since our 
simulations are purely gravitational, with no cosmological 
expansion, the current age of the universe has no direct 
meaning for the simulations themselves. As a result, we use 
the slightly different convergence criterion that the two-body 
relaxation time be longer than the length of time the sim- 
ulation has run. We can invert Eqn. [5] to find r(t r ), the 
minimum radius which satisfies t r > tf for the final time of 
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the simulation tf,we refer to this as our minimum converged 
radius (see Tabic 0. 

Our test runs of single isolated haloes show that the 
density profiles are stable at all radii r(t r ) < r < 0. 87-200 
for the full 5.0 Gyr of simulation: any change in the pro- 
file is under 1%. (Near r2oo there are small deviations due 
to truncating the halo.) Inside r(t r ), the behavior is dic- 
tated by the structure of the halo. NFW haloes have "tem- 
perature inversions" where the maximum velocity disper- 
sion is reached not at the centre but at a finite radius. 
As a result, two-body relaxation does not cause a contrac- 
tion of the central regions as is the case with the singu- 
lar i sothermal sphere (the so-c alled gravothermal catastro- 
phe (IBinnev fc Trem ainc 1987J)) but rather an expansion as 
energy is transfered from regions with higher velocity dis- 
persion to lower. Any model with a cusp shallower than 
7 = 2 starts with a temperature inversion (|Quinlanlll99rt ) 
and therefore undergoes an expansion of the central region 
and a concurrent reduction in central density. The central 
expansion is then followed by a collapsing s tage, initiated 
once the temperatur e inversion is eliminated (|Quinlanl l996: 
lHavashi et al]|2003l) . This process occurs on the evapora- 
tion timescale, which is much longer than the length of our 
simulations at all radii. We do see evolution as a result 
of two-body interactions well within r(t r ) for our isolated 
haloes with 10 5 particles, but seems to be minimal even at 
0.5r(i r ). We conclude that our simulations will be converged 
on all scales greater than r(t r ) and are likely reliable to even 
smaller radii. 

After obtaining stable isolated haloes, we performed 
production runs (all using the massively parallized version 
of GADGET) of halo collisions with different particle numbers 
to check for convergence in the inner structure of merger 
remnants. We did production runs with zero impact param- 
eter using 2 x 10 5 and 2 x 10 6 total particles for the core-core 
collisions, and 10 7 particles (in addition to both lower par- 
ticle number simulations) for an cusp-cusp collision; for a 
core-cusp collision, we did one run with 2 x 10 6 total parti- 
cles. For the two runs with non-zero impact parameter, we 
used 2 x 10 6 total particles. Our two unequal mass merg- 
ers were both core-core collisions with mass ratios of 3 : 1. 
Both runs used the parameters described in Section 2.1 for 
the more massive galaxy. In one run (core4), we reduced all 
characteristic radii by a factor of 3 1//3 , yielding a core radius 
of 5.55 kpc and a cut-off radius of 55.5 kpc. In the other run 
(core5), the less massive galaxy had the same length scales 
as the more massive halo: r CO r C = 8 kpc and r cut = 81 kpc, 
giving a halo that is spatially identical to the heavier one 
but with the density uniformly lower by a factor of 3. By 
using these parameters, we can identify how a merging halo 
having either a smaller core or a lower central density affects 
the merger process. 

For all runs, we use a force softening of 0.175 /i -1 kpc 
for our particles, the new cell-opening accuracy parame- 
ter (errTolForceAcc) of 0.001, and timesteps that are in- 
versely proportional to the acceleration (timestep criterion 
1 of GADGET). (Our choice of errTolForceAcc is somewhat 
smaller than the maximum allowable value suggested by 
IPower et alJ i2003T) in their recent thorough investigation of 
numerical parameters in cosmological simulations, but our 
tests have shown this value is necessary to obtain results 
that are independent of machine architechture and whether 
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Figure 1. Density profiles for 3 different resolution simulations 
of the merger product of two cuspy haloes at 5 Gyr. The profiles 
are indistinguishable at large radii, but the profile for the simula- 
tion with 2 x 10 5 particles (dotted line) shows signs of numerical 
relaxation within ~ 2 /i" 1 kpc. 

or not the run is parallel.) A typical run involved between 
7000 (core-core runs) and 10000 (cusp-cusp runs) timesteps. 
Since the crossing time t c — \J rf^/GM is approximately 
1 Gyr for our haloes, an average timestep corresponds to 
approximately 5.9 x 10 -4 t c . The parameters for all of the 
production runs, as well as some results, are summarized in 
Table 

In all of our simulations, at least 65% of the total ini- 
tial mass ends up within the virial radius of the final halo, 
meaning there are over 6 x 10 6 particles within r2oo in our 
highest resolution simulation, a large enough number to re- 
liably resolve density contrasts exceeding 10 6 . In all runs, at 
least 97% of the particles stays bound to the final halo. All 
mergers appear to be mostly completed by 3.5 Gyr, with the 
virial mass of the haloes increasing by only ~ 10% between 
3.5 to 5.0 Gyr. Fig. shows the density profiles at the fi- 
nal timestep for the three cusp-cusp collisions. The profiles 
of the remnants are extremely similar at all radii greater 
than 2 h~ x kpc, and runs "Cusp2" and "Cusp3" have indis- 
tinguishable density profiles for radii greater than 0.7 h^ 1 
kpc. These results indicate two-body relaxation has not sig- 
nificantly affected any of our haloes at radii larger than 2 
h~ x kpc and that our two highest resolution cusp-cusp colli- 
sions yield density profiles that are reliable on sub-kiloparsec 
scales. 

2.3 Definition of Halo Centres 

The choice of how to define a halo's centre is crucial in an- 
alyzing the final state of the merger remnant since quan- 
tities such as the spherically averaged density profile and 
radial velocity dispersion depend on the location of the cen- 
tre. Previous studies have mainly used two methods for 
defining the centre of a halo: the most bound particle(s) 
(MBP) and the centre of mass (COM). For example, NFW 
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Figure 2. Scatter plot of five snapshots for simulation "Mixed" of a cored and a cuspy halo. The collision is along the horizontal axis, 
and the units of both axes are h~ 1 kpc. Left: both haloes. Middle: cored halo only. Right: cuspy halo only. The halo centres cross at 
~ 1.2 Gyr, then recross at ~ 2.2 Gyr, at which point the relative centre of mass motion is almost completely damped. 
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Figure 3. Time evolution of the kinetic (T), potential (W), and 
total (T + W) energies, as well as the virial ratio (2T/\W\) for 
the cusp-core merger simulation. The virial ratio returns to unity 
soon after the first two encounters at ~ 1.2 and 2.2 Gyr. The total 
energy is conserved to better than 0.4% during the run. The units 
on the energy axis are arbitrary. 

jNavarro. Frenk. fc Whitell996l) use the COM; Moore et al. 
(1999a) use the MBP and note that it gives very similar 
results to the COM; Pearce et al. (1993) use the MBP and 
find that using the COM gives a consistently shallower inner 
profile. 

We adopt an iterative version of the COM technique, 
which we find to be largely independent of initial assump- 
tions. We have tested this method of computing a halo's 
centre against using either the MBP or the COM of the five 
hundred most bound particles. The resulting spherically av- 
eraged density profiles are nearly identical at all radii that 
we consider to be converged - we find no evidence that one 
method systematically yields steeper cusps than the other. 
It is also worth noting that in all cases, the density profiles 
using different criteria for convergence of the centre of mass 
calculation differ only at radii that we already consider to 
be unreliable due to two-body relaxation and Poisson noise. 
As noted earlier, the density profiles are quite similar over 
a factor of 50 in particle number (see Fig. . 



3 RESULTS 

3.1 Scatter Plots; Energy 

Fig-IUshows five snapshots at t = 0, 1, 1.5, 2.5, and 5 Gyr 
from the merger simulation of a cored and a cuspy halo (run 
"Mixed"). (The general sequence of events is the same for all 
three types of mergers listed in Table 1.) To illustrate the 
evolution of the individual haloes, we also plot separately 
the cored halo (middle panels) and cuspy halo (right pan- 
els). During the encounter between 1 and 1.5 Gyr, the outer 
regions of the haloes pass through each other mostly unim- 
peded; the main effect in the outer portions is for the halo to 
become more extended. The panels at 2.5 Gyr show that the 



Table 2. Fits (to Eqn. 1111 for three mergers involving cuspy 
haloes at 5 Gyr. 



Name 


A 


B 


C 




(lO lo A/ kpc- 3 ) 


(h- 1 kpc) 


(cusp slope) 


Cusp2 


0.00593 


9.75 


0.672 


Cusp3 


0.00655 


10.0 


0.696 


Mixed 


0.00425 


12.1 


0.647 



cuspy halo's outer region is more easily stripped away than 
that of the cored profile; this is as expected, since the NFW 
profile is more dense in the centre but less tightly bound 
towards the outer region. 

The merger remnants of the head-on collisions are no- 
ticeably prolate (extended along the collision axis) with axis 
ratios of ~ 1.6:1:1 at half-mass radius. This result is quite 
similar regardless of the initial density profile. The two off- 
set simulations "Core3" and "Cusp4" with orbital angular 
momentum result in mildly triaxial remnants. 

Fig. [3]shows the evolution of the kinetic, potential and 
total energies as well as the virial ratio for the same run in 
Fig. |21 The initial crossing of halo centres is seen to occur 
at ~ 1.2 Gyr, followed by a second crossing at ~ 2.2 Gyr. 
Subsequent oscillations are strongly damped, and the sys- 
tem moves towards a steady state in virial equilibrium. The 
major features in Fig. [3] are common to all production runs 
listed in Table 1. All runs conserve energy to better than 
0.4% and have virial ratios 2T/\W\ within 2% of 1.0 by the 
final timestep. 

3.2 Density Profile 

In Fig. [I] we plot the spherically-averaged density profiles 
p(r) of the merger products at 3.5 Gyr (dashed curve) and 
5.0 Gyr (solid curve) for each of the three types of collisions: 
core-core (left), cusp-cusp (middle), and core-cusp (right). 
The density profiles in all cases are quite stable and constant 
from 3.5 to 5.0 Gyr, indicating that the merging process is 
essentially complete by 3.5 Gyr (also evident in Fig. [3J and 
that our simulations with 2 x 10 6 particles do not suffer 
much from relaxation effects for the radii shown in Fig. 0] 
The dotted curves in Fig. [I] show our fitting results to be 
discussed below. For comparison, we overlay in a separate 
plot (Fig. |5J the final profiles of the merged haloes at 5 Gyr 
from the three collisions. 

The final product of the core-core collision does not 
show signs of moving towards an NFW profile: the remnant 
has a core of approximately the same size as in the origi- 
nal haloes. In order to quantify our results, we have tried to 
fit the final profile to the same functional form as used for 
the initial cored profile (Eqn. [Q. We find the inner region 
of the remnants well fit by a core radius of r coro = 7.8 ft -1 
kpc (for run "Core2"), which is nearly identical to the ini- 
tial r corc = 8 hT 1 kpc, but the density in the outer region of 
the remnants falls off more gradually, as r~ 3 instead of the 
initial exponential cutoff in Eqn. Q. Overall, we find the 
merged halo of core-core collisions better fit by 

p(r) = {r 2 + A B 2 )3/ 2 > ( 10 ) 

with A « 5 x 1O 1O M and B » 11 h' 1 kpc, than by the orig- 
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Figure 4. Density profiles for the core-core (left), cusp-cusp (centre), and cusp-core (right) simulations with 2 X 10 6 total particles, along 
with fits from Eqns , ITC1 and iTTI Both the cusp-cusp and the cusp-core collisions are well fit by an NFW-like profile with an inner cusp of 
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inal profile Q. Since the remnants have virial radii about 
20% larger than their progenitors, this means the relative 
size of the core in the final product is sma ller, a result 
consistent with those of earlier simulations jWhitdll97st 
iPearce. Thomas, fc Couchmanlll993l) . The high spatial and 
mass resolution used here allows us to simulate scales more 
than five times below the core radius and to conclude that 
the core is preserved during one head-on major merger. 

In the cusp-cusp collision, the product is essentially a 
new cuspy halo with an NFW-like form, which we find well 
fit by 



p{r) = A- 



B'' 



(11) 



for radius down to the relaxation scale r = 0.7 kpc for 
the N p = 10 7 particle run. However, the best-fitting value 
for the logarithmic slope of the cusp is shallower than 1, with 
C fa 0.7, perhaps hinting that collisions of galaxies serve to 
reduce the cusp at very small radii. The best-fitting param- 
eters for the density profiles of our two highest resolution 
cusp-cusp collisions, as well as for the core-cusp collision, 
are summarized in Table [5] 

In the case of the core-cusp merger, any hint of the core 
is entirely eliminated. The remnant has an overall density 
profile much closer to NFW than to the original core model 
and again well fit by Eqn. 1111 with an inner slope of C ~ 
0.65. 

To understand the core-cusp merger in more detail, we 
calculate the spherically-averaged density profile for parti- 
cles in each halo separately, first at t — 0.0 Gyr and again at 
t = 5 Gyr; the results are shown in Fig.[S] Interestingly, the 
particles originating from the cuspy halo maintain a cuspy, 
NFW-like profile (left panel), while the particles originating 
from the cored halo retain their original density profile to 
a remarkable extent as well (right panel). The final density 
profile is dominated by the cusp particles in the inner ~ 10 
kpc and by the core particles from ~ 10 — 100 kpc; outside 



10" 



10" 



10" 



10" 



10" 



10 



10" 



-6 




core-core 
. cusp— cusp 
cusp-core 



1 



(h- 



10 

! kpc) 



100 



Figure 5. Density profiles for the core-core (dotted), cusp-cusp 
(solid), and cusp-core (dashed) runs with 2 X 10 6 particles at 5.0 
Gyr. The cusp-core collision results in a remnant with a density 
profile much more similar to the cusp-cusp collision than to the 
core-core one. 



of ~ 200 kpc, both types of particles seem to contribute 
approximately equally (also see Fig. I1UI . 

The results above are all for head-on collisions with a 
zero impact parameter. In Fig. Q we compare the density 
profiles of the non-zero impact parameter runs with the cor- 
responding zero impact parameter simulations. The profiles 
are almost exactly the same, leading us to the conclusion 
that introducing moderate amounts of orbital angular mo- 
mentum into a merging system will not influence our results. 

To ensure that the results are not specific to equal-mass 
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Figure 6. Density profiles of the cuspy particles (left) and core particles (right) initially (think solid) and at 5 Gyr (thick solid) in the 
core-cusp collision simulation. The particles originating from the cuspy (cored) halo clearly maintain the cuspy (cored) profile after the 
merger. The dashed line shows the total density profile of the merged halo at 5 Gyr. 



mergers, we plot the density profiles of our two 3 : 1 mass 
mergers in Fig. |H| along with the density profile of the rem- 
nant in an equal mass merger for reference. The remnants of 
both unequal mass mergers retain cores of ~ 8kpc, the same 
result we found in the equal mass case. In fact, the final den- 
sity profiles in all three are remarkably similar; this shows 
that core-core major mergers produce cored remnants. 

3.3 Structure of Central Region 

Although our simulations are collisionless, resulting in a 
phase space density that must be conserved due to Liou- 
ville's theorem, we have no a priori reason to believe the 
configuration and momentum space distributions will be 
independently conserved, i.e. that a merger of two cored 
(cuspy) haloes will necessarily result in a cored (cuspy) 
remnant. A system's coarse-grained distribution function / 
is constrained to be non- increasing as a function of time 
jBinnev fc Tremaindll987l) . and thus / in the final state is 
required at every point to be less than the maximal initial /, 
but this requirement does not prevent cusps from becoming 
cores or vice-versa in collisionless simulations. In fact, vi- 
olent collisions provide strongly fluctuating potentials that 
cause particles to rapidly exchange energy with the back- 
ground, a process that could quickly affect the phase space 
structure of both initial haloes. Since the main objective of 
our study is to investigate the effects of major mergers on 
the phase-space distributions of dark matter in the inner re- 
gions of galaxy haloes, we examine this region in more detail 
in this section. 

In the core-core and cusp-cusp collisions, we find that as 
expected, both progenitors contribute essentially equally to 
the particles at each radius of the merged halo. As Figs. 2] 
and |H| show, the inner density profile of the merger rem- 
nants is similar to their progenitors, i.e., cored haloes remain 
cored; cuspy haloes remain cuspy. 

The situation is more intriguing for the mixed core-cusp 
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Figure 7. Density profiles for cusp-cusp (thin lines) and core-core 
(thick lines) runs with N p = 2 X 10 6 for b = (solid) and b = 
0.5r200 (dashed) at 3.5 Gyr. The profiles are extremely similar, 
indicating that introducing a moderate amount of orbital angular 
momentum does not significantly impact the merger remnant. 

collision. In order to elucidate the inner structure of this 
type of mergers, we first examine where the initial inner 
particles in each halo end up in the merged halo after 5 
Gyr. (The results at 3.5 Gyr are very similar.) Fig.|J|]shows 
the fraction of particles in each halo that start out within 
radius r initially and also remain within r at 5 Gyr. (By 
construction, this fraction approaches 1 at large r.) It shows 
that in the innermost several kiloparsecs of the merged halo, 
about 60% of the cusp particles (solid curve) that originated 
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Figure 8. Density profiles for three different core-core mergers at 
3.0 Gyr. Dashed line: run Core2 (equal mass, both core radii = 8 
kpc). Dotted line: run Core4 (Mi = 3M2 with core radii n = 8 
kpc, T2 = 5.55 kpc). Dashed line: run Core5 (Mi = 3M2 with 
core radii r% = T2 = 8 kpc). All three mergers produce similar 
remnants with cores that are unreduced relative to the larger of 
the two initial cores, showing that mergers of cored haloes result 
in cored remnants. 
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Figure 10. Ratio of the number of core particles within radius 
r to number of cusp particles within radius r for the cusp-core 
collision at t=0 Gyr (dashed) and t=5 Gyr (solid). The remnant 
is even more centrally dominated by cusp particles than would 
be the case if the two initial haloes were superimposed, further 
showing how the core particles dominate the inner structure of 
the remnant. 
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Figure 9. Fraction of particles starting out within radius r at 
t=0 Gyr that also end up within radius r at t=5 Gyr (from run 
Mixed). The cuspy particles (solid line) are much more likely to 
stay within the inner 10 kpc than the cored particles (dashed 
line). 



from within a given radius stay within that same radius, 
while very few of the core particles (dashed curve) remain. 
We believe it is this "lack of mobility" of the cusp particles 
that gives the cuspy halo a stronger influence at the final 
remnant's centre and forces the final density profile to be 
cuspy rather than cored. 

Next we examine the relative contribution of cusp vs 
core particles to the total density (or potential) in the core- 
cusp collision. Fig. HUI shows the ratio of core to cusp parti- 
cles initially (dashed curve) and at 5.0 Gyr (solid curve). The 
figure shows that the cusp particles dominate the core par- 
ticles number- wise in the centre initially (if the two galaxies 
were superimposed with a common centre), and that at 5 
Gyr, the central region of the merger remnant becomes even 
more dominated by cusp particles. This result implies that 
the tightly bound particles in a cusp tend to stay bound 
together, while a core is much easier to disrupt. This com- 
plements the information from Fig.[5]by showing that for all 
particles of a given type, regardless of origin, the cuspy par- 
ticles tend to take over the central region. Moving outward, 
the core particles begin to contribute more to the remnant. 
Note that the radius where the number of core particles is 
equal to the number of cusp particles stays almost constant 
atrxi 20/i _1 kpc. This is one indication that there is not 
a major reconfiguration of the particles occuring during the 
course of the simulation. We have also calculated the quan- 
tities in Figs. I§1 and llUl at 3.5 Gyr and find them to be very 
similar to those found at 5 Gyr. This corroborates the evi- 
dence in Fig. [3] that the merger is essentially completed by 
3.5 Gyr. 
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Figure 11. Velocity dispersions ay, <re, cr^ (top) and velocity anisotropies (bottom), /3 = 1— (<jg +c^)/2(r^, as a function of halo radius 
r from the cusp-cusp merger run with 10 7 particles (run "Cusp3"). Three time outputs are shown: Gyr (left), 3.5 Gyr (middle), and 
5 Gyr (right). All three velocity dispersions increase after the collision, and the initially isotropic state becomes somewhat anisotropic 
outside the central region. 



3.4 Velocity Profile 

The particle velocities in our initial haloes are isotropic by 
construction, but those in our final haloes all exhibit some 
degree of anisotropy. For the cusp-cusp collision, Fig. 1111 
shows the velocity dispersions (top panels) and velocity 
anisotropy (bottom panels) as a function of halo radius r 
at t = 0, 3.5, and 5 Gyr for the simulation with 10 7 parti- 
cles (run "Cusp3" in Table 1). The velocities remain mostly 
isotropic (i.e. /3 = 0) near the centre of the merger remnant, 
but become mildly anisotropic (f3 < 0.40) with increasing 
radius due to radial infalls. 

We have also examined the velocity distribution func- 
tion for the same cusp-cusp run. We find that in the central 
region of the merged haloes, the velocities are quite Gaus- 
sian in addition to being isotropic; Fig. 1121 shows f(vi) vs. 
Vi for i = r,6,(f> for 10 4 particles in a radial bin arbitrarily 
centred on 10.106 kpc. The distribution function of v r (left 
panel) is well fit by a Gaussian (dashed line) with width 
ay = 177 km s _1 , computed from the particles in that bin. 
The distributions of vg (middle panel) and (right panel) 
are also reasonably well fit by the corresponding Gaussians: 
ag = 153 km s -1 and ct^ = 161 km s _1 . Moving away from 



the centre, we find the radial distribution becomes substan- 
tially less Gaussian and the angular distributions remain 
mostly Gaussian in the intermediate portion of the halo. 
Near the virial radius, however, none of the distributions 
are Gaussian. 

The core-core collision exhibits similar trends as in 
Figs. Ill| and ll2l so we will not show the results here. 

For the core-cusp collision, we plot in Fig. ll3l the radial 
velocity dispersion ay as a function of radius r at t — 0, 2.5, 
and 3.5 Gyr for the core (crosses) and cusp (diamonds) par- 
ticles separately. (Little evolution occurs after 3.5 Gyr.) The 
dispersion of the particles in the cuspy halo changes mod- 
erately between 0.0 and 3.5 Gyr, with the peak in a r rising 
and moving outward but the overall shape maintained. For 
the particles in the cored halo, however, the velocity disper- 
sion near the centre has almost doubled to a r ~ 230 km s _1 
after the merger. 

To understand why the velocity dispersion of the core 
particles increases toward the centre of the merged halo 
upon merging with a cuspy halo, it is instructive to consider 
a situation where the two haloes are spherically symmetric 
and share a common centre. In this case, the potential $ is 
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Figure 12. Velocity distribution functions (solid) for v r (left), vg (middle), and vj, (right) for 10 4 particles in radial bin 10.0371 < r < 
10.1745 kpc from the cusp-cusp run with a total of 10 7 particles at 5 Gyr. The distributions are well-fit by Gaussians (dashed) with 
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Figure 13. Radial velocity dispersion oy as a function of halo radius r from the core-cusp merger run. Three time outputs are shown: 
Gyr (left), 2.5 Gyr (middle), and 3.5 Gyr (right). Particles originating from the cored (crosses) and cuspy haloes (diamonds) are shown 
separately. These two types of particles have comparable a r initially, but the core particles end up with significantly higher oy in the 
inner 10 kpc than the cusp particles after the merger. The cusp particles are also heated, especially at r ~ 10 — 30 h^ 1 kpc, as a result 
of the merger. 



a function of the radial coordinate only and can be written 
as $(r) = $i(r) + $2(j"), a sum of the individual potentials 
$i(r) = $ C orc and $2(r) = "i&cusp. For particles in either the 
cored or cuspy halo, we then have 

it \ 1 f°° /7-.N f d&i d$2\ 

^ = W)L dRp{R) {im + m) 

where M to t(R) ~ Mi(R) + Mi{K) is the total mass interior 



to radius R. From the point of view of the core particles, 
adding the cuspy particles has the effect of adding consider- 
able extra mass and thus deepening the potential. 

If the density profile of the core particles remains 
roughly constant, then adding extra mass in the form of 
cuspy particles essentially heats the core particles to higher 
a velocity dispersion, as seen in Fig. 1131 The reverse is also 
true: the core particles also heat the cusp particles. This 
effect is actually seen in Fig. 1131 as the cuspy particles' dis- 
persion increases relative to the original value, especially 
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Figure 14. Local phase space density, denned as F = p(r)/cr 3 (r), for the N p = 2 X 10 6 simulations of core-core (left), cusp-cusp (centre), 
and cusp-core (right) mergers. The densities are normalized to be equal to that of the cusp-cusp simulation at r = 1 h" 1 kpc. The 
cusp-cusp and cusp-core collisions give phase space densities that can be approximated by a single power law, F oc r~ 1,75 (overplotted 
line), for r < T2Q0 while the core-core collision cannot. 



noticeably near the original peak of ~ 10 h~ kpc. Since the 
cuspy particles dominate the total mass at radii less than 
this, however, the core particles are a small perturbation and 
do not change the cuspy particles' dispersion significantly. 
The change in the peak of the NFW particles' dispersion 
from ~ 10 to ~ 30 h' 1 kpc can also be understood in this 
manner, since the core particles actually contribute more of 
the interior mass at this radius, as shown in Fig. [TU1 

3.5 Phase Space Profiles 

The previous subsections discussed the density and velocity 
profiles separately. It is also instructive to examine the local 
phase space density p(r)/a 3 (r), which is shown in Fig. 1141 
for the three types of mergers studied in this paper. We find 
that for both the cusp-cusp (middle panel) and core-cusp 
(right panel) mergers, although the individual p(r) and a(r) 
are not a pure power-law in r (see Figs. 2] 1111 and 1131 . 
the ratio is interestingly well-approximated by a power law 
p/a 3 oc r~ a , where the best-fitting slope is a « 1.75 for 
both mergers. The core-core merger (left panel), on the other 
hand, shows a core in the phase space density since both p(r) 
and a(r) are constant at r < r covc . 

Our power-law r esults can be compared with 
iTavlor fc Navarrol (1200 ll) , who have suggested that the phase 
space structure of haloes is intimately related to that of 
self-similar collapse onto spherical density perturbations in 
an expanding universe. They solved the Jeans equation un- 
der the assumption of isotropy and power-law p/a 3 oc r~ a 
and found that haloes with the NFW profile are well- 
approximated by the same power law, a = 1.875, as the 
self-similar solutions in secondary spherical infall models 
<Bertschingeilll98Et iFillmore fc Goldreichlll98i lGotdll975t 
lGunnlll977l) . Since the phase space density is inversely pro- 
portional to the local entropy, Taylor & Navarro were led 
to the conclusion that CDM halo structure may be driven 



toward the most sharply peaked phase space distribution al- 
lowable for a monotonically decreasing density profile. For 
a power law phase space density, this solution yields density 
profiles that are very similar to the NFW profile but with a 
slightly shallower inner cusp, p oc r~ ' 75 rather than r _1 . 

Using the power-law result p/a 3 oc r~ 1,75 in Fig. 1141 we 
obtain from the Jeans equation a density profile p oc r -0,72 , 
which is consistent with the best-fitting inner slope shown 
in Sec. 3.2. One can also look at the question in reverse, i.e., 
given a power-law density p oc r -7 , what solutions are ad- 
mitted for the phase space density? Calculations show that 
7 « 0.7 cusps lead to phase space densities well approxi- 
mated by p/a 3 oc r~ 1,75 as r — > 0, but the slope of p/a 3 is 
not strictly a constant beyond the central region (Singh & 
Ma, in preparation). 



4 DISCUSSION 

Our simulations offer insight into the mechanisms that drive 
a merging syste m toward an e quilibrium structure, e.g., vi- 
olent relaxation jLvnden-Belll ll967) and phase mixing. Vi- 
olent relaxation is due to a time-varying potential which 
allows the non-conservation of individual particles' energies. 
The relative motion of the two haloes at the start of our sim- 
ulations provides such a fluctuating potential. From Fig. [3] it 
is clear our mergers involve a rapid transfer of orbital energy 
to internal energy of the individual haloes (and subsequently 
the merger remnants). The decay of orbital motion occurs 
approximately on the crossing timescale, with essentially all 
energy in the form of the remnants' internal energy by 2.5 
Gyr. That the orbital motion is so heavily damped and that 
the remnant haloes do not expand or collapse significantly 
indicate that the time- varying potential necessary for violent 
relaxation is eliminated relatively quickly. One possible ex- 
ception is that density waves could enable violent relaxation 
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Figure 15. Average energy change (in arbitrary units) versus 
impact parameter for the cusp-core collision. Squares are particles 
from the cuspy halo, triangles are particles from the cored halo. 
All particles with small impact parameters tend to lose energy, 
while those with large impact parameters gain energy. The effect 
is more dramatic for cuspy particles than cored ones. 
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Figure 16. Differential energy distribution for the cusp-core 
merger, the solid lines show the cusp particles and the dashed 
lines show the core particles; thin lines indicate the initial distri- 
butions and thick lines the final distributions. Both distributions 
are broadened as a result of the merger, and shell-like structures 
are visible in the final distributions. The energy units are arbi- 
trary. 



to persist even once the haloes have merged. Phase mixing 
tends to damp these waves strongly, however, with a decay 
timescale of r ~ X/a for a perturbation of size A. Taking 
ay = 150 km s _1 to be typical for our merged haloes, this 
means that phase mixing will damp out all perturbations 
smaller than the virial radius within ~ 1 Gyr of the haloes 
merging. In our simulations, this seemingly short timescale 
for violent relaxation is sufficient to drive the velocity dis- 
tribution function in the inner region of the merged halo to 
a Gaussian form , as predicted by statistical reasoning (e.g. 
iNakamuralboOOTi . That the distribution function becomes 
less Gaussian with increasing radius is also in line with ex- 
pectations because much less mixing is possible in the outer 
regions of the halo on the timescales we are studying. 

Following the energy distribution through the collision 
can help in understanding the structure of the final product. 
A particularly interesting quantity is the change in energy 
AE for a particle from initial to final state as a function of 
position. iFunato. Makino. fc Ebisuzakil jl992l) describe an 
energy segregation in their collision simulations: particles 
with higher (less negative) energy tend to gain more en- 
ergy than particles with lower energy. This means a particle 
near the centre of an initial halo should, on average, lose 
energy (i.e. become more tightly bound) relative to a par- 
ticle near the virial radius of a halo. We follow the energy 
change of each particle, and plot in Fig. ll5l the change in en- 
ergy as a function of impact parameter for the core-cusp run 
"Mixed" . As Funato et al. describe, the particles with small 
impact parameters tend to lose energy, while the particles 
with large impact parameters tend to gain it. The effect is 
more dramatic for the cuspy particles than for the core par- 
ticles, perhaps because the inner portion of the cuspy halo 
is initially more tightly bound but wit h an essentially equal 
velocity dispersion profile. As noted by Mcr ral fc Henriksenl 



(2003), however, this energy segregation does not prevent 
the system from equilibrating via violent relaxation, at least 
in the central region. 

Another interesting quantity is the differential energy 
distribution dN/dE, which can be written as a product of 
the phase-spac e distribution function f(E ) and the density 
of states g(E) feinnev fc Tremainelll987l) : 



dN 
IE 



= f{E)g(E). 



(13) 



The shape of dN/ dE is therefore not the same as the shape 
of f(E). In Fie. ll6l we plot dN/dE for the core-cusp merger, 
for both types of particle, at the start of the simulation 
(thin curves) and at 5 Gyr (thick curves). It shows that 
dN/dE is broadened, a natural outcome of mergers where 
the range of energies accessible to an individual particle 
can change drastically. The wave-like features in dN/dE 
at 5 Gyr are indications of shells of pa rticles formed via 
phase- mixing (ISpereel fc Hemauislll992) . We note that un- 
like f(E), which generally decreases with increasing E, the 
shape of dN/dE is skewed towards higher energies because 
of the density of states. Both f(E) an d dN/dE have been 
studied anal ytically for isolated NFW dWidrowll200fj|) and 
King haloes dBinnev fc TremaineHl98^l in equilibrium. The 
overall shape of dN/dE resembles that shown in Fig. 1161 



5 SUMMARY 

We have simulated major mergers of galaxy haloes with both 
cuspy and cored inner density profiles using GADGET, a ver- 
satile N-body tree code. Our findings can be summarized as 
follows: 

(i) Mergers of two cored haloes of equal mass result in 
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a cored halo. The core radius of the merged halo is nearly 
identical to those of the original haloes, but the core to virial 
radius ratio is reduced, with the outer density profile of the 
merger remnant well approximately by p oc r~ i . This result 
also holds for mergers with parent halo mass ratios of 3:1, 
implying all major mergers of two cored haloes result in 
cored haloes. 

(ii) Mergers of a cuspy halo with either a cored halo or a 
second cuspy halo result in a cuspy halo. If the cuspy halo 
starts out with an NFW form, the merged halo also has 
an NFW-like form, with a slightly shallower inner cusp of 
p oc r -0 ' 7 instead of 

(iii) The particles from a given halo retain memory of 
their progenitor's density profile. In major mergers of a cored 
halo with a cuspy halo, the particles originating from the 
cuspy halo end up with a cuspy profile, and those from the 
cored halo end up with a cored profile. (The cuspy particles 
dominate the total potential at the centre of the merged 
halo, so the resulting profile is cuspy.) The core particles in 
the inner 10 hT 1 kpc of the remnant are significantly heated. 

(iv) The merger remnants relax from the inside outwards, 
attaining Gaussian velocity distributions near their centres 
but not throughout the halo. Violent relaxation due to time- 
varying potential is effective only during the initial phase of 
the mergers; phase mixing is likely the dominant relaxation 
process at late times. 

(v) Including a moderate amount of orbital angular mo- 
mentum in the form of a non-zero impact parameter has 
little discernible effect on our findings. 

The set of simulations presented in this paper has been 
performed for mergers of non-spinning haloes in a set bound 
orbit. This clearly only probes a restricted range of param- 
eter space. A possible area for future investigations is to 
include non-zero initial spins, which can influence tidal in- 
teractions significantly but may only affect the final phase 
space structure weakly <IWhitelll979l) . We have also chosen 
not to convolute the issues in this paper by ignoring minor 
mergers with satellite haloes and continuous accretion of dif- 
fuse material. The latter mechanisms may be responsible for 
generating the first generation of cuspy haloes from initially 
cored haloes, a process not achieved by major mergers of 
cored haloes according to our study. 

Although the length and time scales quoted in this pa- 
per are for Milky- Way sized haloes of 1O 12 M0, there is noth- 
ing in the numerical code or in our analysis that is special to 
this particular mass scale (except a mildly mass-dependent 
concentration parameter for NFW haloes). Our conclusions 
in this paper are thus likely to be applicable to major merg- 
ers of comparable-mass haloes over a wide mass range. Our 
finding that cuspy haloes are resilient to major mergers hints 
that the existence of cuspy haloes is natural in the hierar- 
chical structure formation scenario: the build-up to larger 
mass haloes via repeated merging may be sufficient to pro- 
duce cuspy haloes from cored haloes and preserve the general 
structure of cuspy haloes. 
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